Meshless Chebyshev RPIM Solution for Free Vibration of Rotating Cross-Ply Laminated Combined Cylindrical-Conical Shells in Thermal Environment

This paper provides a numerical solution to the vibration of a rotating cross-ply laminated combined conical-cylindrical shell in the thermal environment. Its numerical discrete solution method uses the meshless method. The combined shell assumed the temperature independence of material property is divided to the fundamental conical and cylindrical shell substructures, and the theoretical formulation for each substructure is derived based on the first order shear deformation theory (FSDT) and Hamilton’s principle. The effects of the initial hoop tension and temperature change are considered through the kinetic energy reflecting the effects of centrifugal and Coriolis forces and additional strain energy by the nonlinear part of the Green–Lagrange strains. The substructures are then assembled according to the continuity conditions. The boundary and continuity conditions are simulated by introducing artificial virtual spring technology. The displacement component in the theoretical formulation is approximated using a meshless Chebyshev-RPIM shape function. The reliability of the method is verified by comparing with mature and reliable results. The free vibration characteristics of the rotating combined conical-cylindrical shell structure under various sizes, speeds and temperatures are given by numerical examples.


Introduction
In the aerospace field, laminated shell structures are widely used in the shell structures of gas turbines and high-power aircraft engines [1][2][3][4][5]. In these high-end fields, the vibration of the structure will bring huge economic losses, so it is necessary to study its free vibration behavior before designing such a structure.
With the progress of computational science, many different methods such as the Haar wavelet discretization method [6], geometric analysis (IGA) method [7], spectral-Tchebychev solution technique [8], Ritz method [9,10] and finite element method [11][12][13] are employed for dynamic characteristics analysis of the composite structures. Ye et al. [14] derived the classical open shell formula on the basis of FSDT, and used Chebyshev polynomial to construct the displacement shape function, and solved the free vibration frequency of the open shell through the Rayleigh Ritz program. Caresta and Kessissoglou [15] reported a wave solution for the free vibrational frequencies of a homogeneous composite conical-cylindrical shell, where the displacement component was approximated by a power series. Tornabene et al. [16] reported a method for dynamic analysis of laminated hyperbolic shells and rotating panels on elastic foundations using GDQ. Li et al. [17] reported a Jacobi Ritz method for solving the free vibrations of laminated hyperbolic rotating shells with general boundary constraints. In addition, in recent decades, several studies on the dynamic mechanical properties of rotating structures in thermal environments have been developed. Shakouri et al. [18] reported the vibrational behavior of a conical shell of a functionally graded material with temperature-dependent material properties during rotation. Afshari [19] extended the generalized differential quadrature method to the solution of free vibration of a rotating conical shell reinforced by graphene nanomaterials. Bhangale et al. [20] reported a finite element method for analyzing the dynamics of functionally graded conical shells operating in high temperature environments. Tian et al. [21] obtained the free vibration and forced vibration solutions of the combined conical cylindrical shell by the dynamic stiffness method. Qin et al. [22] used the Rayleigh-Ritz method which based on the energy variation principle to solve the free vibration problem of a cylindrical shell-ring-plate coupling system. Singha et al. [23] analyzed the free vibration characteristics of rotating pretwisted sandwich conical shells in a thermal environment based on high-order shear deformation theory by using a finite element method. Talebitooti et al. [24] investigated the frequency behaviors of the joined conical-conical panel structures based on FSDT by applying Hamilton's principle. Soureshjani et al. [25] investigated the free vibration behaviors of composite joined conical-conical shell in the thermal environment by using a generalized differential quadrature method. Shi et al. [26] proposed an analytical model for investigating the vibration characteristics of a functionally graded conical-cylindrical coupled shell structure by using a spectro-geometric method. Ghasemi et al. [27] investigated the influences of distribution, mass and volume fractions of fiber, boundary conditions and lay-ups on the sensitivity of vibration behaviors of hybrid laminates cylindrical shell according to Kirchhoff Love's first approximation shell theory. Liu et al. [28] focused on the influences of rotation on the frequencies and critical speed of CNTs/fiber/polymer/metal laminates cylindrical shell based on Love's first approximation shell theory. Semnani et al. [29,30] analyzed the vibration behaviors of microshell under varied working conditions by using a finite element method.
In addition to the above methods, the development of meshless theory provides a brand-new idea for plate-shell vibration analysis. Based on the three-dimensional elastic theory, Kwak et al. [31,32] proposed a meshless strong-form solution for the free vibration of laminate shells. In their method, Chebyshev polynomials are introduced as basis functions in the construction of shape functions. In the meshless approach, the establishment of the system algebraic equations of the problem domain does not use a pre-defined mesh for domain discretization, but instead uses nodes [33]. Zarei et al. [34] constructed the displacement function of a prestressed laminate using meshless radial basis point interpolation and analyzed its vibration characteristics. Mellouli et al. [35] used the same method to build a vibrational analysis model of functionally graded carbon nanotube-reinforced shells. Zhang et al. [36] introduced the vibrational behavior of carbon nanotube-enhanced functionally graded triangular plates using a meshless method. Fallah and Delzendeh [37] studied the free vibration of laminates with meshless finite volume method (MFV) as the model solution method and moving least squares approximation to approximate the displacement component. Kwak et al. [38] combined the Chebyshev polynomial with the radial basis point interpolation method to construct the displacement shape function of the open laminated cylindrical shell with elliptical section, and solved its natural frequency.
The purpose of this paper is to study the vibration properties of a rotating crosslaminated conical-cylindrical shell in a thermal environment using meshless theory, considering that the combined structure is divided into cylindrical shell and conical shell structure, and the cylindrical shell is a special conical shell. Therefore, the equations of motion suitable for rotating conical shells are first established within the FSDT framework. Then, the two substructures are assembled by the continuity equation to obtain the equation of the overall structure. The effects of centrifugal force, Coriolis force, and temperature are considered in the equations of motion, and the displacement components involved are approximated using a meshless TRPIM shape function. The accuracy and reliability of the proposed method are verified through the convergence study and comparing with the results of the literature and ABAQUS. Finally, the effects of parameters such as geometry, temperature difference and rotational speed on the free vibration of the cross-ply laminated composite conical-cylindrical shell structure are studied. To sum up, the investigations of this paper can analyze the variation tendency of vibration characteristics of rotating crosslaminated conical-cylindrical shell in the thermal environment and provide the theoretical basis for the designation and manufacture of rotating cross-laminated conical-cylindrical shell structures which are used in aircraft, missiles, submarines, etc. Figure 1 shows a model of laminated combined conical-cylindrical shell rotating with rotating angular velocity Ω under the influence of temperature difference ∆T. The symbols L 1 and L 2 denote the lengths of the two meridians. The thickness of the combined shell is uniformly set to h. ϕ represent the semi-vertex angle of conical shell. The symbols R 1 and R 2 represent the radii at both ends of the conical shell, respectively. The cylindrical shell is connected at the big end of the conical shell, so the radius of the cylindrical shell is also R 2 . The orthogonal curvilinear coordinate system (x, θ, z) is introduced into the middle surface of each substructure. The orthogonal coordinate system o-xθz is established on the middle surface of the substructure, then the radius R of the random position on the conical shell is as follows:

Governing Equations and Boundary Conditions
According to the assumption of first-order shear deformation [39], the displacement (u, v, w) of any position on the elastic structure can be represented by the displacement (u, v, w, ψ x , ψ θ ) of the mid-plane.
Combined with the linear elasticity theory, the relationship between the stress and displacement of the shell is defined: represents the normal strain and shear strain of the elastic element, and χ = χ α , χ β , χ αβ T represents the bending and torsional curvature changes of the elastic body. γ = γ 0 βz , γ 0 αz T denotes transverse shear strain. A and B denote the Lamé parameters.
The matrix form of the stress resultants-strain relationship of moderately thick crossply conical shell is as  where N = N α , N β , N αβ T , the internal element represents the in-plane force. M = M α , M β , M αβ T , the element represents the bending moment, and Q = Q β , Q α T is the shear force vector. A represents the tensile stiffness matrix, B is the bending stiffness matrix, and D is the coupled tensile bending stiffness matrix. A c denotes the shear stiffness matrix. Their specific form is: where N denotes the number of laying layers of the laminate, k c = 5/6 is the shear correction coefficient and the symbol Q k ij denotes the elastic stiffness coefficient [38]. In the thermal environment, the thermal stress of the kth layer in the cross-layer is expressed as follows: where the thermal expansion coefficients α k ij of the kth layer are given by where δ k and α ij denote the fiber angle of the kth layer and linear thermal expansion coefficients along the principal axes of a layer, respectively. The thermal strain can be written as nonlinear part of Green-Lagrange strain.
Substituting Equations (2) and (14) into Equation (13) ε In the thermal field, the strain energy of the structure is expressed as: Meanwhile, when the shell rotates, an initial hoop tension will be generated due to centrifugal force, which will generate a part of the strain energy.
where N 0 θ = ρhΩ 2 R 2 is the initial hoop tension, the unit of ρ is kg/m 3 . The kinetic energy is When the shell is not affected by external force, according to the variational principle, the equilibrium equation and boundary conditions of the heated rotating cross-layer shell are deduced.
The obtained governing equations are expressed as: where the matrices C and M are expressed as: where the inertia terms are The boundary conditions obtained from Hamilton's principle are expressed as: Combining Equations (21), (24) and (25), we can derive the governing equations and boundary conditions for the cross-layer shell in the thermal physics field.
where ω and n denote the natural frequency and circumferential wave number, respectively.
Substituting Equation (20) into Equations (21) and (24), the one-dimensional governing equations and boundary conditions of rotating cross-ply conical shell in thermal environment are obtained.

Meshfree TRPIM Shape Function
The radial point interpolation method is a newly developed meshless method, which is an important and widely used method for solving partial differential equations. The unknown displacement function u(x) is approximated by using the RPIM difference of the polynomials and can be defined as in [38].
where R i (x) is the radial basis function (RBFS), and n r is the number of nodes of the point x in the support domain. p j (x) is the polynomial in the space coordinate x T = (x, y), and n p represents the number of polynomials. If n p = 0, it is a single radial basis function (RBFS), otherwise it is an RBF with n p polynomial basis functions added. Generally, for a onedimensional problem, the basis function of the polynomial is p j (x) = [1,x,...,x np ] T , and in a two-dimensional problem, the polynomial basis is p j (x) = [1,x,y,...,x np ,xy np−1 ,...,yx np−1 ,y np ] T . However, using a power function polynomial basis is often inaccurate in solving differential equations. Chebyshev polynomials have important applications in approximation theory. Corresponding interpolation polynomials minimize the Longo phenomenon and provide the best consistent approximation of polynomials in continuous functions. Therefore, this study uses Chebyshev polynomials as interpolation basis functions. where The multi-quadrics (MQ) radial function with shape parameters αc and q are used in this paper.
where r i denotes the distance between the supported point x J (J = 1,2,n r ) in the supported domain and calculated node x I . For the one-dimensional problem in this paper r i = |x Jx I |, d c is a characteristic length related to the node spacing in the support domain of the compute node. When the nodes are evenly distributed, d c is the distance between adjacent nodes. Otherwise, dc is the average node spacing within the node distribution domain.
In meshless theory, the size of the local support domain will affect the interpolation accuracy, and a suitable size of the supported domain should be selected [32]. The size of the supported domain of the calculated node can be characterized as follows.
where α s represents the scale factor of the support domain.
In order to determine the coefficient vectors a and b in Equation (30), a support field for calculated node x I needs to be formed, which includes n r field nodes. Let Equation (30) satisfy the calculation of n node values around point x I , which yields n r linear equations. The matrix of these equations can be expressed as the following form.
where R 0 represents the RBFs matrix and T nt is the Chebyshev polynomial matrix [30]. The coefficient vector a of RBFs is expressed as follow. a = a 1 a 2 · · · a n r T The coefficient vector b of the Chebyshev polynomial basis function is written as follow: Since there are n r + n t unknowns in Equation (35), a unique solution cannot be obtained, so it is necessary to add n r equations through the following constraints to make the coefficient matrix of the equation system full rank.
Combining Equations (35) and (38), the matrix representation of the following system of equations can be generated.
Delete unnecessary terms in the Chebyshev-RPIM shape function above, and obtain the Chebyshev-RPIM shape function corresponding to the final node displacement.
Through the above derivation, the displacement components of the nodes can be expressed as follows.

Discretization of Governing Equations and Boundary Conditions
The substructure of the combined structure is discretized using N nodes, and the displacement approximation function at node x I is represented by a Chebyshev-RPIM shape function.
where N s represents the number of nodes covered by the support domain, I 5 represents a 5 × 5 identity matrix. Substituting Equation (47) into Equation (26) to obtain the discretized governing equation represented by node information.
where the nodal matrices K xI , C xI and m I are as follows.
Similarly, the discrete equations for whole system are obtained by assembling those of each node according to the node number [40].

Continuous Condition
The governing equations and boundary equations of the substructure have been deduced and discretized before, but a complete solution system has not been established. The combined structure can be divided into conical shell and cylindrical shell. According to their geometric characteristics, considering their displacement continuity and physical coordination, the right boundary of the conical shell and the left boundary of the cylindrical shell can be modified as follows.
where k b denotes the connection stiffness between substructures, and symbols co and cy denote conical and cylindrical shells, respectively. The matrix form of the continuous condition can be written as follows.
B xco Φ T co U sco + K 12 Φ T cy U scy = 0 : Right boundary of conical shell B xcy Φ T cy U scy + K 12 Φ T co U sco = 0 : Left boundary of cylindrical shell (50) where U sco and U scy are the displacement vectors of the nodes of the cylindrical shell and the conical shell on the coupling interface, respectively. The coupled stiffness matrices K 12 and K 21 are as follows.
Finally, matrix assembly is performed to obtain the vibration control equation of the overall structure. ( The natural frequency of the conical-cylindrical composite structure in the thermal environment is obtained by the harmonic response method.

Numerical Results and Discussions
This paper provides a meshless free vibration analysis model of a rotating combined conical-cylindrical shell structure in a thermal environment. The proposed method is compiled with MATLAB software. The number of nodes and the size of the support domain will affect the convergence effect of the algorithm. After obtaining the appropriate support domain size and number of nodes through convergence analysis, the numerical results are compared with finite element software or published literature to ensure the reliability and accuracy of the proposed method. Then, focusing on structural characteristic parameters and the effect of external physics on structural frequencies, some parametric study cases are provided. Unless otherwise stated, the natural frequencies of the considered combined shells are expressed in the dimensionless parameters as ω * = ωR 1 ρ/E 2 and the material properties of the layers are given as: E 1 = 175 GPa, E 2 = 32 GPa, µ = 0.25, G 12 = G 13 = 12 GPa, G 23 = 5.7 GPa, ρ = 1760 kg/m 3 , α 11 = 1.2 × 10 −6 , α 22 = 2.3 × 10 −6 and α 12 = 0. The symbols C, F, and S are used to represent the tightened boundary conditions, free boundary conditions and simply supported boundary conditions, respectively. The corresponding boundaries are described as follows: C: k = k v = k w = k x = k θ = 10 14 . S: k u = k v = k w = k θ = 10 14 . k x = 0, F: k u = k v = k w = k x = k θ = 0. Then define boundary rules. For example, CF represents that the boundary of the conical shell segment is a fixed boundary, and the boundary of the cylindrical shell segment is a free boundary.

Verification and Convergence Study
First, according to the basic theory of the meshless method, the key factor affecting the convergence of numerical results is the number of nodes. Therefore, before the numerical comparison and parametric analysis, the advanced convergence analysis is carried out to ensure that the obtained calculation results are stable. Table 1   In the previous convergence analysis, it has been determined that the meshless theory is applied to the structural vibration analysis, and the obtained results have good stability.
However, it has not been demonstrated whether the obtained results have a high level of confidence. Therefore, it is necessary to further compare the results obtained in this study with the existing publications or the results obtained by the finite element software ABAQUS. Table 2 compares the vibration frequencies of the non-rotating combined conicalcylindrical shell, considering no temperature difference between the inside and outside of the shell. The dimensions of the structure are:R 1 = 0.4226 m, ϕ = π/6, L 2 = R 2 = 1 m and h = 0.01 m. The material properties are: E = 211 GPa, ρ = 7800 kg/m 3 , µ = 0.3. The dimensionless frequency of a non-rotating combined conical-cylindrical shell structure is defined as: ω * = ωR 2 ρ(1 − µ 2 )/E. The results obtained by the meshless method are compared with the published literature [10] and [15], and the difference between the results obtained by the meshless method and the literature is very small. Table 3 compares the frequency results obtained by different numerical methods for rotating isotropic combined conical-cylindrical shells. The boundary conditions, geometry and Poisson's ratio of the combined structure are the same as those in Table 2, and the rotational speeds considered are 0.01 rad/s, 100 rad/s, and 500 rad/s, respectively. The comparison results show that the method in this paper is in good agreement with the results in the literature. Finally, it is verified that the model established in this paper can be applied to the structural vibration solution in a thermal environment. In Table 4, the vibration frequency of the nonrotating laminated combined conical-cylindrical shell structure in a thermal environment is analyzed using the finite element software ABAQUS and the method in this paper, respectively. The considered structural geometry is: The temperature change is 50K. The frequencies obtained by these two methods are in good agreement. Figures 2 and 3 represent the mode shapes of laminated combined conical-cylindrical shell, corresponding to the natural frequencies from Table 4. Meanwhile, it is necessary to point out the fact that the following numerical discussion illustrates that this method can be used to analyze structural vibration behavior in thermal environments. All in all, after sufficient comparison, it is proved that the method established in this paper can be applied to the vibration analysis of the rotating composite conical shell and cylindrical shell in a thermal environment.

Numerical Examples
In Section 3.2, we conduct sufficient comparative verifications to demonstrate that this method can analyze the vibrational behavior of the rotating composite cone-column structure in a thermal environment. First, the effect of semi-vertex of the combined conicalcylindrical shell structure with T = 0 K on the natural frequency is studied. The structure shape is: It can be seen from Figure 4a,c,d that as the half-apex angle of the conical shell increases, the frequency of the combined structure gradually increases slightly first and then decreases significantly under the CC, SS and CS boundary conditions, respectively. At the same time, the difference between the forward wave frequency and the backward wave frequency of the conicalcylindrical composite shell with different rotational speeds is getting smaller and smaller, and the influence of rotational speed is also weakened. As on can see from Figure 4b, under the CF boundary condition, the backward wave frequency of the composite structure corresponding to Ω = 150 rad/s and Ω = 200 rad/s will decrease first and then increase with the increase of the half apex angle, and the rest of the natural frequency change curves all decrease. Likewise, with the same rotational speed, the gap between the forward traveling wave and the backward traveling wave of the structure also decreases.
Secondly, Figure 5 studies the variation of forward wave frequency with temperature for a combined conical-cylindrical shell with a rotational speed of 50 rad/s. In Figure 5a,c,d we selected the forward wave frequency with the circumferential wave number n = 1~4 and the axial half-wave number m = 1 as the research object to discuss its variation with ∆T. The variation interval of the temperature difference is [0 K, 500 K], and the boundary conditions of the studied structures are CC, CF, SS and CS boundary conditions, respectively. The structural geometry parameters are the same as in Figure 4. It is clear from Figure 5b that for the combined structure under the CC, CF, SS and CS boundary conditions, respectively, when n = 0 and n = 1, the curve representing the relationship between frequency and temperature difference approaches the horizontal line. At this time, the effect of temperature difference on them is minimal. However, for n = 2~4, the natural frequency of the structure decreases with increasing temperature difference. Under the CF boundary condition, for n = 2 and n = 3 with the increase of the temperature difference, the frequency value first decreases, and then does not change. At this time, at the turning point of the curve, the structure undergoes thermal buckling.  In addition, the effect of the rotational speed on the frequency of the combined structure in the thermal environment is also studied, where three temperature differences are selected as 0 K, 200 K and 1000 K, respectively, the rotational speed variation interval is 0,500], and the boundary condition is CC, CF, SS and CS, respectively. Following the structure of Figure 5 for analysis, the research results are shown in Figure 6. Meanwhile, it is necessary to point out that the ω * f (n=1,m=1) and ω * f (n=3,m=1) are considered in the following discussion. As can be seen in Figure 6a,c,d, the frequency of all intercepted backward wave frequencies under CC, SS, CS boundary conditions, respectively, increases as the rotational speed increases. However, under CC, CF, SS, CS boundary conditions, respectively, the forward wave frequencies for ω * f (n=1,m=1) decrease and for ω * f (n=3,m=1) increase with increased rotating speed. As shown in Figure 6b, under CF boundary conditions, the variation tendencies of forward and backward wave frequency are the same as the above other boundary conditions. However, the forward and backward wave frequencies produce model jumping with temperature differences increased, and the above phenomenon can be attributed to thermal buckling.
Finally, Tables 5 and 6 show the dimensionless frequencies of rotating cross-ply combined conical-cylindrical shell with various geometry and boundary condition in thermal environment. The geometrical parameters of the structure and the temperature difference are given in the table header, and the lamination scheme is [0 • /90 • /0 • ]. It can be seen from Table 5 that as the length of the cylindrical shell increases, the stiffness of the structure decreases, and both the forward wave frequency and the backward wave frequency gradually decrease. The rules in Table 6 are the same as those in Figure 4, and thus are not repeated here. These results are valuable to designers and serve as benchmarks for future numerical studies.

Conclusions
This paper focuses on the free vibration analysis of laminated combined conicalcylindrical shell rotating in a thermal environment. The equations of motion of the whole system are derived by combining the equations of individual substructures obtained using Hamilton's principle, in which the effects of temperature change, centrifugal and Coriolis forces are taken into account. For numerical calculation of equations of motions, the meshfree strong form method using TRPIM shape function is applied. Through the convergence study, the number of node is determined. The accuracy and reliability of the proposed method are verified through comparison with the results of finite element program. Finally, the effects of parameters such as geometric dimensions, rotating speed and temperature change on the free vibration of combined conical-cylindrical shell are investigated through some numerical examples. The conclusions obtained in this study are as follows: (1) The meshless Chebyshev-PRIM technique is effective and has relatively high accuracy in the vibration solution of rotating structures. This method has the advantage of fast convergence, and relatively accurate results can be obtained with a smaller number of nodes. (2) The increase of the half-apex angle of the conical shell reduces the structural rigidity, so the structural frequency decreases. For the combined structure under the CC boundary, after the cone angle increases to a certain extent, the effect of the rotational speed will decrease, and the frequencies corresponding to different rotational speeds will gradually approach. (3) If the temperature is too high, thermal stress is accumulated inside the structure, the stiffness of the structure is reduced, and the frequency of the combined structure will also decrease. For the boundary conditions with weakened constraints, such as the CF boundary, thermal buckling also occurs with the increase of the temperature difference.